Vortices and Ring Solitons in Bose-Einstein Condensates 



O 

o 

(N 

o 
O 



L. D. Carri'2 and Charles W. Clark^ 
1. Physics Department, Colorado School of Mines, Golden, CO 8O4OI 
2. Electron and Optical Physics Division, National Institute of Standards and Technology, 
Technology Administration, U.S. Department of Commerce, Gaithersburg, Maryland 20899 

(Dated: February 2, 2008) 

The form and stability properties of axisymmetric and spherically symmetric stationary states in 
two and three dimensions, respectively, are elucidated for Bose-Einstein condensates. These states 
include the ground state, central vortices, and radial excitations of both. The latter are called ring 
solitons in two dimensions and spherical shells in three. The nonlinear Schrodinger equation is taken 
as the fundamental model; both extended and harmonically trapped condensates are considered. It 
is found that instability times of ring solitons can be long compared to experimental time scales, 
making them effectively stable over the lifetime of an experiment. 
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I. INTRODUCTION 

One of the primary motivations in the original deriva- 
tion of the Gross-Pitaevskii equation was to describe vor- 
tices in superfluids P, 0, |1] . The quantization of vortic- 
ity is a central difference between classical fluids and su- 
perfluids 0. The Gross-Pitaevskii equation has proven 
to be an excellent model for weakly interacting, dilute 
atomic |l|i[l|l|i,[l3 and molecular |il|Tl[ia Bose- 
Einstein condensates (EEC's), which are generally su- 
perfluid. In the context of optics the Gross-Pitaevskii 
equation is known as the nonlinear Schrodinger equa- 
tion (NLSE) [13. Vortices Il5l Il6l , as well as their one- 
dimensional analog, solitons |l7l Il8l | , have been observed 
in EEC's in numerous experiments. The NLSE describes 
these observations well [lll20ll2lll2^. 

However, there are in fact much richer vortex and soli- 
ton structures yet to be observed. One such structure 
is the ring soliton a soliton 

extended into two dimensions which loops back on itself 
to form a ring, i.e., a radial node. This appears as an ax- 
isymmetric radial node of the condensate; set at the right 
distance from the origin, it becomes a stationary state. 
In the optics context, the ring soliton is well known to be 
unstable to the snake instability, whereby it decays into 
vortex-anti-vortex pairs. In this article, we not only de- 
scribe vortices and radially excited states of EEC's with 
high precision, but we also show that in a harmonic trap 
the decay time can be long compared to the classical os- 
cillation period of the trap and even the lifetime of the 
condensate itself |30l |. 

We consider three cases for the external, or trapping 
potential V{r). First, V{r) = corresponds to an in- 
finitely extended condensate and leads to solutions of 
beautiful mathematical form. Second, V{f) = V{r) = 
for r < R and V{r) = 00 for r > R, corresponds to a disk 
in two dimensions and a sphere in three, i.e., an infinite 
well. This introduces confinement into the problem, con- 
nects heuristically with the first case and general knowl- 
edge of solutions to Schrodinger equations, and serves as 
a bridge to the experimental case of a harmonic trap. 



Third, F(r) = ^M[uj^{x'^ + y^) + iv^^z^], with w < w^, 
where M is the atomic mass and lij,uJz are the classi- 
cal oscillation frequencies, corresponds to a highly oblate 
harmonic trap, which is most relevant to present experi- 
ments. 



This work follows in the spirit of a previous set of in- 
vestigations of the one- dimensional NLSE, for both re- 
pulsive and attractive nonlinearity |3ll l32| . In that work 
it was possible to obtain all stationary solutions in closed 
analytic form. In the present cases of two and three di- 
mensions, we are unaware of an exhaustive class of closed 
form solutions but instead use a combination of analyti- 
cal and numerical techniques to elicit the stationary and 
stability properties of similar solutions. Here, we treat 
the case of repulsive atomic interactions; as in our previ- 
ous work on one dimension, the attractive case has been 
treated separately [s^ , due to the very different character 
of the solutions. 



Equations similar to the NLSE are often used as mod- 
els for classical and quantum systems. Thus a tremen- 
dous amount of theoretical work has been done on vor- 
tices, to which the reader is referred to Fetter and 
Svidzinsky on EEC's Donnelly on Hehum II [3, and 
Saffman on classical vortices [s^ as good starting points 
for investigations of the literature. As NLSE-type equa- 
tions apply in many physical contexts, our results are 
widely applicable beyond the EEC. 



The article is outlined as follows. The derivation of the 
fundamental differential equations is presented in Sec. ITU 
In Sec. mil the ground state and vortices in two dimen- 
sions are presented. In Sec. llVl the stationary radial ex- 
citations of these solutions arc illustrated. In Sec. Ivlthe 
ground state and its radial excitations in three dimen- 
sions are treated. In Sec. IVII the stability properties of 
solution types containing ring solitons are discussed. Fi- 
nally, in Sec. IVIII we discuss the results and conclude. 
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II. FUNDAMENTAL EQUATION 



Constant Potential 



The fundamental differential equation is derived as 
follows. The NLSE, which models the mean field of a 
BEC 0,1113, is written as 



d 



where V{r) is an external potential, g = A-Kh^as/M, 
Us is the s-wave scattering length for binary interac- 
tion between atoms with Us > 0, since this is the re- 
pulsive case, and M is the atomic mass. The conden- 
sate order parameter ip = ip(r,t) = yjn{r^ t) exp[iS'(r, t)], 
where n(r, t) is the local atomic number density and 
v{f,t) = {h/My^ S{r,t) is the local supcrfluid velocity. 
Note that, in two dimensions, the coupling constant g is 
renormalized by a transverse length [H, Ho, 113, Ell • 

We assume cylindrical or spherical symmetry of both 
the external potential and the order parameter in two 
or three dimensions, respectively. This has the effect 
of reducing Eq. to one non-trivial spatial variable. 
Specifically, wavefunctions of the form 

V'(^) = fmir) exp{im4)) exp{-i^d/h) exp{i6a) (2) 

are treated, where m is the winding number, /i is the 
eigenvalue, also called the chemical potential, (j) is the 
azimuthal coordinate, r is the radial coordinate in two or 
three dimensions, and 0o is a constant phase which may 
be taken to be zero without loss of generality. 

Assuming an axisymmetric stationary state in two di- 
mensions of the form given in Eq. Eq. becomes 



2M \ dr'^ r dr / 
+gfl + V{r)U - f^fm = . 



(3) 



Assuming spherical symmetry in three dimensions, one 
finds a similar equation. 

In Eq. it is assumed that to = 0, in keeping with the 
spherical symmetry. This is an important special case, 
as it includes the ground state. In the remainder of this 
work, TO will be taken as non-negative since = f\m\- 

The physically relevant solutions of Eqs. ||2Jl and ^ in- 
clude the ground state, vortices, ring solitons, and spher- 
ical shells, as we will show in the following three sec- 
tions. A variety of solution methods are used to treat 
the wavefunction in different regions of r, as discussed in 
Appendix fXl Different rescalings of Eqs. ©-lEl) are ap- 
propriate for the three potentials we consider: constant, 
infinite well, and oblate harmonic. These are treated 
bricfiy in the following subsections. 



A constant potential has no direct experimental real- 
ization. However, it does have mathematical properties 
which are helpful in understanding the confined case, not 
to mention beautiful in themselves. For instance, the vor- 
tex solution manifests as a boundary between divergent 
and non-divergent solutions, as shall be explained. 

The potential V{r) = Vq can be taken to be zero with- 
out loss of generality. Then the variables can be rescaled 
as 

^ (5) 
X ^ • (6) 




Note that the radial coordinate is scaled to the length 
associated with the chemical potential. Then Eq. Q 
becomes 



1 



and Eq. Q becomes 



II , ' 3 I A 

X 



(7) 



(8) 



where /i has been assumed to be positive, since this is 
the physically meaningful case for repulsive atomic inter- 
actions. Note that, in these units, the length scale of a 
vortex core is on the order of unity. 



B. Infinite Potential Well in Two and Three 
Dimensions 

Let us first consider the two dimensional infinite well. 
A confined condensate may be obtained by placing an in- 
finite potential wall at fixed r — -\- y'^ , at any node of 
the wavefunction ip{r, (p). This treats either a condensate 
tightly confined in z or a cylinder of infinite z extent. In 
either case, one derives a 2D NLSE from the 3D one by 
projecting the z degree of freedom onto the ground state 
and integrating over it. This leads to a straightforward 
renormalization of the coefficient of the nonlinear, cubic 
term. Extremely high potential wells in the x-y plane 
have been created in BEC experiments via higher order 
Gauss-Lagucrre modes of optical traps [s^l . 

For the purposes of our numerical algorithm outlined 
in Appendix it is convenient to keep the same scalings 
as Eqs. (jSJ-ljni. Then the normalization has to be treated 
with some care. The normalization of ip in two dimen- 
sions for a BEC of N atoms in an infinite cylindrical well 
of radius R is given by 



drr\il}{r,(t);t)f . 



(9) 
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After the change of units given by Eqs. (jSJ and ^ the 
normalization becomes 



where 



M 

M= :^gN, 



(10) 

(11) 
(12) 



are the effective nonhnearity and cyhnder radius. 

The three-dimensional infinite well is treated similarly, 
with the normalization being J\f, rather than 2M, due to 
the extra angular integration. 

The details of the algorithm for calculation of quan- 
tized modes in two and three dimensions is discussed in 
Appendix IX4l 



C. Oblate Harmonic Trap 

For the harmonic potential we shall focus on the case 
of a highly oblate trap, which is the experimentally rel- 
evant one to obtain axial symmetry in two effective di- 
mensions |4C| . Then, again projecting onto the ground 
state in z, the rescaled nonlinearity Af has a simple in- 
terpretation: 



TV = 2asNy/Mu}z/2Trn , 



(13) 



where lj^ is the angular frequency of the trap in the z 
direction. The trap is isotropic in the remaining two 
directions. All energies can then be scaled 

to hoj, lengths to ^ = y/h/Muj, etc., as follows: 



1 f+if _!^f 

r, \ Jm' ~Jm -2 

2 V r 



where the tildes refer to harmonic oscillator scalings. Ex- 
plicitly, f = r/£, fm = £fm, jj- = ii/htu, and the normal- 
ization is 



drr\f„ 



--Af. 



(15) 



The main difference between Eqs. (|14ll and Eq. is 
that an extra parameter must be set in the numerical al- 
gorithm of Appendix fXl i.e., the rescaled chemical poten- 
tial. Then the normalization is determined from Eq. (|15|l . 




FIG. 1: (color online) Approaching the vortex solution. 
Shown is the dependence of the wavefunction on the precision 
in the critical determining coefficient a'^^^'' for a quantum 
vortex stationary state of the nonlinear Schrodinger equation 
in two dimensions, (a) The vortex solution forms a boundary 
between convergent (ao — fio^'(l — 10^ ^''), solid black curve) 
and divergent (ao = aQ^^(l + 10^^*"), dot-dashed blue curve) 
solutions. As the precision is reduced, the first node moves 
towards the origin and the solution approaches the Bessel 
function: (b) 8 digits of precision; (c) 2 digits of precision. In 
(b) and (c) the regular Bessel function solution to the linear 
Schrodinger equation is shown for comparison IIH (dashed 
red curve). Note that all axes are dimensionless. 



III. GROUND STATE AND VORTICES 

In order to solve Eqs. (0) and we use a numer- 
ical shooting method, as discussed in App. lA 31 Two 
initial conditions arc required, as the NLSE is second or- 
der. These are provided by a Taylor expansion around 



X = 0, as described in Sec. lA II By this method a single 
parameter is sufficient to determine the solution. This 
parameter, ag, is the lowest non-zero coefficient in the 
Taylor series. High precision is required, as discussed in 
the appendix, with the number of digits of precision de- 
termining where the numerical method either converges 
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or diverges. 

For the first potential we consider, ^(x) = 0, the 
ground state, which is obtained for m = 0, hes precisely 
on the border between convergence and divergence of the 
algorithm. The value of ao which is exactly on this bor- 



(a) 



(b) 



der we term 



,(0) 



The wavefunction 770 (x) is constant 



and, according to our scalings, is simply ±1 

Then a^^^ = ±1 and the precision is infinite. Values of 
|ao| larger than 1 lead to a divergent solution, while val- 
ues of |ao| less than 1 lead to a convergent one. In the 
latter case the wavefunction oscillates an infinite num- 
ber of times and approaches zero, as will be discussed in 
Sec. IIVI Thus, in general, I770I can approach only three 
asymptotic values: 0, 1, and 00. 

For the case of non-zero winding number, one finds 
a central vortex. As in the case of the ground state, 
lim^^^oo |?7mP = n = fi/g is the asymptotic density of 
the vortex state. In this region the spatial derivatives 
yield zero and 77m(x) ±1 as x ~* oo- The vortex again 
lies on the border between divergence and convergence 
of our algorithm, given by a single parameter Oq™'' which 
determines the whole Taylor series. In Appendix 1X1 the 
precision issues are discussed in detail. In Fig. ^ is illus- 
trated the algorithmic approach to the vortex solution for 
zero external potential. The effect of the number of dig- 
its of precision is shown in detail, and further explained 
in Appendix^ where the best values of Cq™'' for winding 
numbers m = 1 to m = 5 are given. 

For the second potential, an infinite well in two dimen- 
sions, the properly normalized ground state and vortex 
states are produced by the simple prescription give in Ap- 
pendix ^31 An example is shown for m = 1 in Fig. |2Ia). 
The ground state for three dimensions can be found by 
a similar method, and is shown for m = in Fig.EJa). 

For the third potential, a strongly oblate harmonic 
trap, the form and stability of both the ground state and 
vortex solutions in a harmonic trap have already been 
thoroughly studied elsewhere [I^, IH, 113 ■ Our algorithm 
reproduces all the relevant results of these authors from 
the non-interacting to the Thomas-Fermi regime, as we 
explicitly verified in detail in the case of Ref. ■ 



IV. RING SOLITONS 



Ring solitons can be placed concentrically to form a 
stationary state. In an extended system a denumer- 
ably infinite number are required, as already indicated 
in Fig. n and discussed in further detail in Sec. IIV Cl be- 
low. For an infinite well in two dimensions, radially quan- 
tized modes are distinguished by the number of concen- 
tric ring solitons, as discussed in Sec. IIV Al and illustrated 
in Fig.[5| Note that, for attractive nonlinearity, the num- 
ber of rings can vary from one to infinity, even for an 
infinitely extended condensate, as discussed in Ref. [s^. 
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FIG. 2: Ring solitons in an infinite well. Shown is the form 
of the wavefunction 771 for (a) the ground state, (b)-(e) the 
first four axisymmetric excited states, and (f ) a highly excited 
state, all as a function of the radial coordinate r scaled to 
the cylinder size R. Cylindrical box boundary conditions in 
cylindrical polar coordinates are assumed. The central vortex 
has winding number m = 1. The quantized modes (b)-(f) 
correspond to increasing numbers of concentric ring solitons. 
Here, the case of strong nonlinearity is illustrated, with 27V = 
402. Note that all axes are dimensionless. 



Quantized Modes in the Cylindrical Infinite 
Well 



In Fig. 121 are shown the wavefunctions for the ground 
state and first three excited states for a fixed strong non- 
linearity J\f in the presence of a singly quantized central 
vortex, i.e., winding number m ~ I. The weaker the 
nonlinearity and the larger the number of nodes in the 
solution, the closer it resembles the regular Bessel func- 
tion Jmix)- In appropriately scaled finite system, 
the relative weight of the kinetic term to the mean field 
term in the NLSE increases strongly with the number 
of nodes. We note the contrast with the distribution of 
nodes for the corresponding ID NLSE [sil l33 |. where 
the nodes are evenly spaced even in the extremely non- 
linear limit. The value of the nonlinearity was chosen to 
be 2JV = 402. For transverse harmonic confinement of 
angular frequency oj^ = 2tt x 100 Hz and ®^Rb, which 
has a scattering length of Qs = 5.77 nm, this corresponds 
to ~ 10, 000 atoms. Note that, since rj is scaled to 
g/ and ^ depends on the number of nodes, the vertical 
scaling is different in each panel of Fig. [3 
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FIG. 3: (color online) Ring solitons in an infinite well. Shown 
are the eigenvalue spectra fin as a function of the normaliza- 
tion 2J\f. The winding numbers (a) m = 0, (b) m = 1, and (c) 
m = 2 are illustrated for the ground state (solid black curve), 
and the first three excited states (blue dashed, green dot- 
dashed, and red double-dot-dashed curves). The circles show 
the actual data points. Note that all axes are dimensionless 
as indicated in the text, and on the same log- log scale. 



The eigenvalue spectra for the ground state and the 
first two excited states are shown in Fig. with /i scaled 
to the radius of the infinite well, 



Mi? ^ 



(16) 



The winding numbers m = 0, m = 1, and m = 2 are illus- 
trated on a log-log scale. Clearly there are two regimes. 
For small J\f, hr is independent of the norm. This must 
be the case near the linear-Schrodinger-cquation regime, 
since must approach the eigenvalues of the regular 



Bessel function Jm{x) which solves Eq. lO with no cubic 
term. One finds {^ir)j where x^j""""" is the 

known value of the j'^ node of the Bessel function 
and j also refers to the j^^ quantized mode. 

For large Af, the figure shows that fXR ex JV. This 
dependence can be understood analytically in the case 
of the Thomas-Fcrmi-like profile for the lowest energy 
state with winding number m = 1. Consider the scaling 



p = r/R 
V i V 



(17) 
(18) 



Then Eq. © becomes 



1 dn 
pdp 



-K-inR^UK^ + pnn = 0. (19) 



The Thomas-Fermi profile is obtained by dropping the 
derivatives: 



KTF (p) 



AttRW 



where 



Pm 



(20) 
(21) 

The normal- 
(22) 



Then the chemical potential in units of the energy asso- 
ciated with the cylinder radius is 



is the core size and ktf is zero for p < p„i 
ization condition is 

27ri?2 /■ dpp[KTF{p)f = 1. 



PR 



4AA 



1 - Pm + Pm In(Pm) ' 



(23) 



The limit pm <C 1 is consistent with the Thomas-Fcrmi- 
likc profile, which neglects the radial kinetic energy. In 
this limit, one finds 



W [1 + Pm-Pm In (p^O] . (24) 



For < 0.03, as is the case in the right hand side of 
Fig. |2Ia)-(c) where p^ ~ 10"^, the dependence on is a 
less than 1% perturbation. 



B. Quantized Modes in the Strongly Oblate 
Harmonic Trap 

In figure 21 is shown the wavefunction for fixed strong 
nonlinearity for the ground state and first three excited 
states. As in Sec. lIV Al adding nodes to the wavefunction 
drives the system towards the linear regime, where the 
solution is Bcssel-function-likc. In Fig. |31 arc shown the 



6 




FIG. 4: Ring solitons in a harmonic trap. Shown is the form 
of the wavefunction /i for (a) the ground state and (b)-(d) 
the first three axisymmetric excited states, all as a function 
of the radial coordinate f. The harmonic trap is strongly 
oblate, so that it is effectively two-dimensional. The central 
vortex has winding number m = 1. The quantized modes 
(b)-(d) correspond to increasing numbers of concentric ring 
solitons. Here, the case of strong nonlinearity is illustrated 
for (a), with M = 100; for three rings, i.e., (d), the solution 
already appears nearly linear. Note that the tildes signify 
that all axes are in harmonic oscillator units. 



spectra for the ground state and first three excited states, 
for a winding number of m = 0, 1, 2. Note that the chemi- 
cal potential is rcscalcd to the harmonic oscillator energy. 
As in Sec. IIV Al there are two regimes. For large non- 
linearity a Thomas-Fermi approximation can be applied 
to obtain the asymptotic dependence of the chemical po- 
tential on the nonlinearity, similar to the procedure of 
SecHvH 
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C. Asymptotic Behavior for Zero External 
Potential 

For X ^ oo and |ao| < Ioq"'^! the wavefunction ap- 
proaches zero in an infinitely extended system. Thus one 
expects that the nonlinear term 77^ in Eq. |(7J) is negligi- 
ble in comparison to the other terms, and the differential 
equation returns to the usual defining equation for the 
Bessel functions. The asymptotic form of the regular 
Bessel function Jm{x) is 113 



X 



mix 
~2~ 



(25) 



to leading order in the amplitude and the phase. How- 
ever, one cannot neglect the effect of the cubic term on 
the phase shift, as may be seen by the following consid- 
erations. 

The asymptotic form of the Bessel function can be de- 
rived via the semiclassical WKB approximation 01 ■ The 



FIG. 5: Ring solitons in a harmonic trap. Shown are the 
eigenvalue spectra as a function of the normalization M, 
all in harmonic oscillator units. The ground state (solid black 
curve) and the first three isotropic excited states (blue dashed, 
green dot-dashed, and red double-dot-dashed curves) are il- 
lustrated for a winding number of (a) m = 0, (b) m = 1 and 
(c) m — 2. The circles show the actual data points. Note 
that all axes are dimensionless and on a log-log scale. 



phase shift of 7r/4 can be derived by analytical continua- 
tion or other means |45| . The semiclassical requirement 
that the de Broglie wavelength be small compared to the 
length scale of the change in potential is not quite sat- 
isfied near the origin. The rescaling y = yo ln(x) suf- 
fices to map the problem onto the usual semiclassical 
one. One can avoid the rescaling by the substitution 
rri^ m? — j. Then the term m7r/2 in the phase shift 
follows directly 23| ■ We use this simpler method in 
order to derive the phase shift in the nonlinear problem. 
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The semiclassical momentum is 



where the effective potential is 

2 



(26) 



(27) 



Taking the nonhncar term as pcrturbativc. to lowest or- 
der Eq. (P7|l becomes 



Kff(x) 



B2 



H cos^ (6'„,) , 

X" X 

TTITT TT 



(28) 
(29) 



where i? is a constant coefBcient of the amplitude of 
the wavefunction. In the linear case, it is convention- 
ally taken as B = i/2/7r. Expanding Eq. (^f)]! for large 
X, one finds 



p{x) - 1 - 



cos2 e.„ 



2,Y 

The semiclassical form of the wavefunction is l44l 14 



^ cos - - 



to leading order, where 



5* = 



dx'p(x') 



(31) 



(32) 



is the semiclassical action. Upon substitution of Eq. H3U|) 
one finds the form of the wavefunction to leading order 
in the amplitude and the phase. 



B 

I 



Hx) 



-— - - -f (5(ao,TO) 



(33) 



where 5 is a phase shift which depends on the deter- 
mining coefficient and the winding number m. This 
coefficient cannot be analytically determined by Eq. H32() 
since the large x form of the wavefunction was used, while 
the phase shift is due to its behavior in the small x re- 
gion. The amplitude coefficient i? is a free parameter, the 
square of which is related to the mean number density. 

The form of Eq. H33(l resembles that of the Coulomb 
function in that it has a ln(x) dependence in the 
phase. This is due to the 1/x term in the effective po- 
tential in Eq. H28|l . It is in this sense that the nonlinear 
term in Eq. ^ cannot be neglected, even as 77 0. 



V. SPHERICAL SHELLS 

Spherical shells are the three-dimensional analog of 
ring solitons. Quantized modes in the confined case in- 
volve successive numbers of nested nodal spherical shells. 
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FIG. 6: Spherical shell solutions. Shown is the form of the 
wavefunction for (a) the ground state, (b)-(c) the first two 
excited isotropic states, and (d) a highly excited isotropic 
state. The boundary conditions are an infinite spherical well 
in three dimensions. The winding number is zero. The case 
of strong nonlinearity is illustrated, with J\f = 251. Note that 
all axes are dimensionless. 



In this section, we will treat the two cases of zero poten- 
tial and an infinite spherical well. A power series solution 
of Eq. JSJ may be developed by substitution of Eq. |IA1|I . 
This leads to a solution similar to that of Sec. lA II All 
coefficients in the power series are given as polynomi- 
als in the determining coefficient ag, which is the value 
of the wavefunction at the origin. The special solution 
%(x) ^ 0,0 = il is the ground state in an extended 
system. Thus a^^^ = ±1. Positive values of which 
are larger than unity lead to a divergent solution. Those 
less than unity lead to a convergent solution which ap- 
proaches zero as X — ^ oo- 



A. Quantized Modes in the Spherical Infinite Well 

Solutions can be quantized in the three-dimensional 
spherical well in the same way as the two-dimensional 
case. The solution methods are identical to those of 
Sec. IIV AI In Fig.Elare shown the ground state, the first 
and second excited states, and a highly excited state. A 
fixed nonlinearity oi M = 251 was chosen. In Fig. [7| 
is shown the eigenvalue spectra on a log-log scale. As in 
the two-dimensional case of Fig.|31 there are two regimes. 
For small A/", the eigenvalues arc independent of the non- 
linearity, {^iR)j (x"""^')^- The constant x""''" is the 
distance to the jth nodes of the spherical Bessel function 
jo(x) which solves Eq. ||SJ) with no cubic term For 
large M, one again finds a linear dependence. A simple 
estimate based on the Thomas-Fermi profile for m = 0, 
which is just fo{r) = [ijg for r < R and zero other- 
wise, gives the chemical potential of the ground state as 
Mi? = AA/3. 
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FIG. 7: Spherical shell solutions. Shown are the eigenvalue 
spectra = {2M /TL^)fi as a function of the normahzation 
A/". The ground state (solid black curve) and the first three 
isotropic excited states (blue dashed, green dot-dashed, and 
red double-dot-dashed curves) are illustrated for a winding 
number of m = 0. The circles show the actual data points. 
Note that all axes are dimensionless and on a log- log scale. 



B. Asymptotic Behavior 

As r — > oo the spherical shell wavefunction ap- 
proaches zero in an infinitely extended system. Just as 
in Sec. llV Cl one can use the WKB semiclassical approx- 
imation method to determine the asymptotic form of the 
wavefunction. The solution to Eq. | |H| w ithout the cubic 
term is the spherical Bessel function [4l| 



Jo(x) 



X 



(34) 



where we have assumed the wavefunction to be finite at 
the origin. One can take the nonlinear term as perturba- 
tive since, for sufficiently large the linear form of the 
wavefunction must dominate. Then the effective poten- 
tial in the WKB formalism is 



Vesix) - B — — 



(35) 



where B is the amplitude of the wavefunction. From 
Eqs. (123 and the WKB action is 



S ~x 



4 X 



(36) 



for large where the phase shift (5'(ao, m) cannot be de- 
termined from the large x behavior of the wavefunction. 
Then the asymptotic form of the wavefunction is 



Vo[X) — sm 
X 



(37) 
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FIG. 8: (color online) Stability Properties of Ring Solitons. 
Shown is the Boguliubov linear excitation frequency as a func- 
tion of the nonlinearity A/", all in harmonic oscillator units. 
The index q signifies the winding number of the Boguliubov 
mode, (a) A ring soliton with no central vortex; for TV 25 
the primary instability mode is g = 2. (b) A ring soliton in 
the presence of a central vortex of winding number m = 1; 
for A/" ^ 100 the dominant instability is g = 3. The solu- 
tions are always formally unstable, though instability times 
can be much longer than experimental timescales for small 
nonlinearities. The circles show the actual data points. 



VI. STABILITY PROPERTIES 

The stability properties of ring solitons in a strongly 
oblate harmonic trap are of particular importance, 
as such solutions may be realized in present experi- 
ments [23. We perform linear stability analysis via the 
well-known Boguliubov de Genncs equations [20I l48l| : 



where 



2M 



W^ + V{r) + 2g\if-t,, 



(38) 
(39) 



(40) 



and is the eigenvalue. In Eqs. (|38|I - H39|I 
a complete set of coupled quasiparticlc amplitudes that 



Uj and Vj are 
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obey the normalization condition 

Jd^r{\u,f~\v,f)^l, (41) 

with D the number of effective dimensions. These ampH- 
tudes represent excitations orthogonal to the condensate 
ip. They have a straightforward quantum mechanical in- 
terpretation in terms of a canonical transformation of the 
second-quantized Hamiltonian for binary interactions via 
a contact potential of strength g ■ Quasiparticles are 
superpositions of particles (creation operators) and holes 
(annihilation operators). Classically, they can be inter- 
preted simply as linear perturbations to the condensate. 

If the eigenvalue flj is real, the solution ip is stable. 
If Qj has an imaginary part, then "0 is unstable. There 
arc significant subtleties in Boguliubov analysis; see the 
appendix of Ref. for an excellent discussion. In the 
present effectively 2D potential with a condensate solu- 
tion of the form given in Eq. (01, it is useful to redefine 
the Boguliubov amplitudes as suggested by Svidzinsky 
and Fetter IH: 



u{r) \ 

v{r) J e 



Vq{r) 



(42) 



where f = y^z^^ + y^/f, C is the harmonic oscillator length 
discussed in Sec. Ill CI and we will neglect perturbations 
in the z direction, due to the strongly oblate trap. Equa- 
tion (|42|l represents a partial wave of angular momentum 
q relative to the condensate. Then, in harmonic oscillator 
units (see Sec lIICjl . Eqs. (I2Hll-(|SSll become 



C^Va 



\fn 



I/. 



— ? 

n. 



(43) 
(44) 



where 



1 / 

2 \d?^ 
+2\fm\'' - 



Id (q± mY 9 
— f 



(45) 



The different centrifugal barriers inherent in C± show 
that the two amplitudes behave differently near the ori- 
gin, with Uq oc ■pl'?+™l and Vq cx f'*^™' as f ^ 0. Note 
that the nonlinear coefficient is absorbed into the nor- 
malization of fm - see Eq. (I15|) . 

The condensate wavefunction can be obtained via the 
shooting and Taylor expansion methods described in Ap- 
pendix ^ Then Eqs. (g^-lg^ can be solved straight- 
forwardly with standard numerical methods. We use a 
Laguerre discrete variable representation |5fl l5ll | , which 
is particularly efficient for this geometry, allowing us to 
go to hundreds of basis functions. The winding number q 
of the Boguliubov modes was checked for g = to g = 6 
over the entire domain of our study. 

The results are shown in Fig. |S1 In panel (a), it 
is apparent that a single ring soliton placed on top of 



the ground state, i.e., with no central vortex, is always 
formally unstable to the quadrupole, or g = 2 mode. 
The instability time is given by T = — u;/Im(f2) in har- 
monic oscillator units. Typical trap frequencies range 
from = 27r X 10 Hz to 27r x 100 Hz. Therefore, when 
|Im(17)/ti;| ^ 1, T can be much longer than experimental 
time scales of 100 ms to 1 s. In this case, we say that the 
solution is experimentally stable. The instability time T 
can even be longer than the lifetime of the condensate, 
the latter of which can range from 1 to 100 seconds. Ac- 
cording to Fig. |SJa), this occurs for small nonlinearities. 
For larger nonlinearities, other modes, such as the oc- 
topole {q = 3), also become unstable. In Fig. |S{b) is 
shown the case of single ring soliton in the presence of a 
singly quantized central vortex, i.e., m = 1. The solution 
is again formally unstable, though first to octopolc rather 
than quadropole perturbations. For small nonlinearities 
it is experimentally stable. 

We did not quantitatively study solution stability for 
an infinite well. However, we expect that the boundary 
provides additional stability of ring solitons. To decay, 
ring solitons must break up into pairs of vortices via a 
transverse oscillation. To oscillate, the ring soliton has 
to move away from the barrier and inwards towards the 
origin. This requires shrinking the circumference of the 
ring, which costs energy, as a single ring feels an effective 
potential which pushes it outwards, as for example in an 
unbounded system. We expect that very long decay times 
follow. We make the conjecture that, in the infinitely 
extended system, the presence of an infinite number of 
rings, tightly pressed up against each other, has the same 
effect with regards to the inner ring as the boundary in 
the confined system. 



VII. DISCUSSION AND CONCLUSIONS 

Soliton trains in one-dimensional BEC's, which are 
similar to the nested ring solitons which form radial 
nodes, have been found to be archetypes of planar soli- 
ton motions encountered in three-dimensional BEC's. 
These solutions to the one-dimensional NLSE are sta- 
tionary states with islands of constant phase between 
equally spaced nodes. When ap prop riately perturbed, 
they give rise to soliton motion |3ll IH^ . In fact, the 
stationary solutions can be considered to be dark soli- 
tons in the limiting case of zero soliton velocity, and the 
perturbation that produces propagating "gray" solitons 
is the imposition of a slight phase shift in the wavefunc- 
tion across a node. These one-dimensional examples were 
found to have experimentally accessible analogs in three- 
dimensional BEC's in which optically induced phase 
shifts across aplane of symmetry resulted in soliton mo- 
tion 0, 0, [53 ■ Correspondingly, the two-dimensional 
ring solutions presented herein suggest the possibility of 
creating ring soliton motion by imposing a phase shift 
across the boundary of a disk. We will address this in 
subsequent work. 
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The existence of ring dark soHtons has been predicted 
theoretically [2^ |23| and demonstrated experi- 

mentally [23, in the context of nonlinear optics. It 
has been suggested that a single ring dark soliton could 
be created in a confined BEC ■ A ring dark soliton 
corresponds to a single node in our ring solutions. It is 
known that a single ring dark soliton in an infinitely ex- 
tended system expands indefinitely j2J|. This therefore 
clarifies why the ring solutions require an infinite num- 
ber of nodes in order to obtain a stationary state. It 
also explains why, for cylindrical box boundary condi- 
tions, the creation of nodes tends to be towards the edge 
of the condensate. In Ref. [23 it was found that a sin- 
gle ring dark soliton was unstable to vortex pair creation 
via the transverse, or snake instability in the near-field 
(r <^ -Rtf) Thomas- Fermi approximation in a harmonic 
trap. In agreement with this work, we have shown, with- 
out approximations, that linear instability times can be 
made so long that a ring soliton is in fact stable over the 
lifetime of the experiment. This result holds independent 
of the trap frequency in the 2D plane. 

The ring solutions that we have discussed in Sec. IIVI 
might be realized in an experiment that approximates a 
deep cylindrical potential well, e.g. an optical trap using 
a blue-detuned doughnut mode of a laser field |54l |. In 
such a system, the ground state vortex solution will re- 
semble that of Fig.[3Ia), where the abscissa is the radial 
coordinate in units of the well radius. The first radi- 
ally excited vortex solution will then resemble that of 
Fig. |2Ib) . By use of an optical phase-shifting technique 
such as that employed in Refs. [i3j[i3' might be able 
to generate this state and observe its subsequent motion. 
We showed that the same qualitative pattern of radial 
nodes that occurs for the infinite well is also found in 
strongly oblate harmonic traps. 

Concerning the central vortex core of the ring solu- 
tions, we note that single vortices are quite long lived 
compared to experimental time scales 0, [s^ [S^ • It 
is possible that forced excitation of the condensate may 
couple resonantly, either directly or parametrically, to 
ring formation. The same possibility exists for spherical 
shell solutions in the observation of nodal spherical shells. 
In two dimensions, unlike in three, multiply charged vor- 
tices do not dynamically decay into singly charged vor- 
tices with the addition of white noise to the system [s^ . 
despite their being thermodynamically unstable. In fact, 
Pu et al. showed that stability regions recur for large 
nonlinearity in two dimensions for to = 2 Recent 
experiments have been able to create and manipulate vor- 
tices of winding number greater than unity in a variety 
of ways [5^ [53, 113 ■ Therefore our study of vortices in 
two dimensions of winding number higher than unity is 
experimentally relevant, despite their being thermody- 
namically unstable [63 . 

In summary, we have elicited the form and properties 
of stationary quantum vortices in Bose-Einstein conden- 
sates. It was shown that their axisymmctric stationary 
excitations take the form of nodal rings, called ring soli- 



tons. Quantization of these states can be attained in con- 
fined geometries in two dimensions; we considered both 
the infinite well and a harmonic trap. Similar methods 
were used to study the ground state and isotropic sta- 
tionary excitations in a spherical infinite well. Two im- 
portant aspects of these solutions is that (a) the rings or 
spherical shells pile up near the edge of the condensate, 
rather than being evenly spaced in r, in contrast to the 
one-dimensional case, and (b) the chemical potential de- 
pends linearly on the atomic interaction strength when 
the mean field energy dominates over the kinetic energy, 
i.e., in the Thomas-Fermi limit. We showed that the ring 
solitons are experimentally stable for weak nonlinearity. 

This work was done in the same spirit as our previ- 
ous articles on the one-dimensional nonlinear Schrodinger 
equation |3ll . l33 | . A companion work [33 treats the at- 
tractive case, which has features radically different from 
the present study. For instance, vortex solutions arc not 
monotonic in r. Moreover, there are a denumerably in- 
finite number of critical values of the determining coef- 
ficient ao for fixed winding number which correspond to 
the successive formation of nodes at r ~ 00^ even without 
an external trapping potential. 

Finally, we note that phenomena similar to the spher- 
ical shell solutions have been experimentally observed in 
BEC's lU, while 2D BEC's appropriate to the investi- 
gation of ring solitons are pre sently under intensive in- 
vestigation in experiments [6^ [13 . 
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for useful discussions. LDC thanks the National Sci- 
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of Basic Energy Sciences via the Chemical Sciences, Geo- 
sciences and Biosciences Division for support. The work 
of CWC was partially supported by the Office of Naval 
Research and by the National Science Foundation. 

APPENDIX A: NUMERICAL METHODS AND 
PRECISION ISSUES 

1. Analytic Structure of the Solutions 

The following numerical methods are discussed explic- 
itly for the infinitely extended condensate, i.e., for a con- 
stant external potential. A slight modification for the 
infinite well is discussed in App. lA 41 and briefiy for the 
strongly oblate harmonic potential in Sec. Ill CI 

Since Eqs. {Tj) and © do not contain any non- 
polynomial terms, one may begin with a power series 
solution by a Taylor expansion around x = of the form 

00 

where the aj are coefficients. For solutions which have 
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the limiting behavior ry™ — > at the origin, which 
is necessarily true for all non-divergent solutions with 
TO ^ 0; the nonlinear term rj^ becomes negligible as 
X — > 0. Then the Bessel function solutions to the lin- 
ear Schrodinger equation are recovered. This will equally 
be true where rj has a node, in the neighborhood of the 
node. Thus near the origin the wavefunction must be- 
have as rjmix) x™- This motivates the choice of the 
exponent of x in Eq- I|A1|I . By examination of Eqs. ^ 
and l|SJl, it is clear that only even or odd powers of x 
can have nonzero coefficients. The Taylor series has been 
written in such a way as to eliminate all terms which 
are obviously zero. Substituting Eq. (|A1|) into Eqs. (O 
and JHl, the coefficients can then be obtained recursively 
by equation of coefficients of equal powers. One finds that 
all coefficients aj for j ^ can be expressed as a poly- 
nomial in odd powers of ao of order 2([j/(to + 1)J) + 1, 
where \jx\ denotes the greatest integer less than or equal 
to X. For example, the first few terms for to = 1 are 



ai 

0,2 
04 



192 



9216 
1 

737280 



(ao + 8al) , 
(ao + SOal) , 



(ao -f- 656al + 1152ag) . 



(A2) 



Thus the coefficient oq is the only free parameter of the 
problem. We consider only ao > 0, since for each solution 
?7,ri(x; ao), there is a degenerate solution 7?m(x; ~«o)- 

The power series provides a useful, practical method 
for propagating the solution of the NLSE away from the 
singular point at r = [64| . However, it is not a practi- 
cal method for extension to large r and we therefore use 
other methods in intermediate and large r regions. An 
asymptotic expansion which is not formally convergent 
but nevertheless useful is obtained by the transformation 



Then Eq. {Tj) becomes 



d_ 



2/-2 

m Q 



1 



(A3) 



77™(C)=0. (A4) 



A Taylor expansion around C — yields the asymptotic 
power series solution 



7/1 = 1- 

772 = 1 - 
f?3 = 1 - 
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24661 


2^ 






16x6 


128 x* 
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6 


68 


1514 










x' 


" ■ ■ ■ ; 


9 


153 




4473 


962037 


2^ 






16 X*^ 


128 x« 



,(A5) 
(A6) 
(A7) 



etc. as X 00. Since this series has no free parameters, 
it is clear that only one value of the determining coeffi- 
cient ao in Eq. IjAip can lead to the vortex solution. We 



define this critical value as . All values of ao > a^Q^^ 

lead to divergent solutions, while values of ao < Uq"^ lead 
to solutions which asymptotically approach zero, as shall 
be discussed in Sec. IIVI It is in this sense that the vortex 
solution manifests as a boundary between divergent and 
non-divergent solutions. 

Lastly, it is worthwhile to mention a limiting case 
which is useful in Sec IIVI For to = all coefficien ts a, 
are zero except for ao. Examination of Eqs. O and (|A4|I 
shows that the solution must be 



^(x) 



1 , 



(A8) 



This is the ground state in an extended system in two 
dimensions. 



2. Solution by Pade Approximant 

It is desirable to determine the behavior of the vortex 
in intermediate regions between zero and infinity. The 
two point Pade approximant is defined by the rational 
function 



•A/"p,g(x) _ Co + ClX + h Cq-lX' 



q-l 



where 



do + dix + • • • + dqx'^ 



oixn 



^pAx)~fix)-DpAx) 

as X ^ a-nd 

KAx) ^ 9{x)t^pAx) = oix^-''^-') 



(A9) 



(AlO) 



(All) 



as X ^ 00 for all p such that < p < 2q, with p and q 
integers. The functions /(x)i dix) are power series ex- 
pansions of the same fimction as x 0, 00. The solution 
of Eqs. l|ITO|) and (|ITT|l for the power series expansions 
of Sec. lA ll leads to a determination of the critical value of 
the determining coefficient ap™'* and therefore a solution 
of the NLSE valid over all space. 
For instance, taking (7 = 3, one finds 



m 



V2x + 2x2 
1 + \/2 + 2x2 



(A12) 



The approximation can be successively improved. Taking 
q = 4 one finds 



16V62 X + 248x2 + 30V62 x 



124 + 31 V62 X + 248x2 + -f 30 V62 chi^ 



(A13) 



and so on. By solving Eqs. IjAlOp and (|Alip at suc- 
cessively higher order, one obtains a convergent value of 

(m) 

However, in practice this procedure is limited in pre- 
cision due to the appearance of spurious roots as well 
as by computation time. Due to the nonlinear nature 
of Eqs. l|ITO)l and l|ITT|l . the re are multiple values of ao 
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TABLE L Results of 2-point Pade approximant for the 
determining coefficient ag""' of quantum vortex stationary 
states. flQ™' is the first nonzero coefHcient in the power series 
solution of the 2D NLSE under the assumption of a single 
central vortex of winding number m, and determines all 
subsequent coefficients. 



Winding number m 


(m) 


Precision 





1 


oo 


1 


0.583189 


6 


2 


0.15309 


5 


3 


0.02618 


4 


4 


0.00333 


3 


5 


0.0002 


1 



TABLE II: The best possible converged values via numerical 
shooting methods of the determining coefficient for quantum 
vortex stationary states. 



Winding number m 


"o 


Precision 





1 


oo 


1 


0.5831 8949 5860 3292 7968 


20 


2 


0.1530 9910 2859 54 


14 


3 


0.0261 8342 07 


9 


4 


0.0033 2717 34 


8 


5 


0.0003 3659 39 


7 



which satisfy them. These roots become sufficiently close 
to each other so as to mislead root-finding algorithms. 
Roots close to Og™^ tend to produce solutions which are 
asymptotically correct but have spurious non-monotonic 
behavior, i.e., wiggles in intermediate regions. There are 
two options for root finding of large systems of coupled 
polynomial equations. One may find all roots and test 
them one by one. However, the computation time be- 
comes prohibitive for higher order polynomials and large 
numbers of simultaneous equations. Or, one may use 
Newton's method or some other local root-finding algo- 
rithm to find the root closest to the correct one for the 
previous lowest order. The appearance of spurious roots 
then becomes a limiting factor. 

In Tabic m is shown the best convergent values of Oq"'' 
for winding number zero to five. Higher m leads to the 
appearance of more spurious roots and therefore a lower 
maximum precision. This may be understood as follows. 
The coefficients in the power scries defined by Eq. (|AH) 
were polynomials in ao of order 2(\_j / [m -f- 1)J ) -f 1. High 
winding number therefore requires a greatly increased 
number of terms in order to obtain improved values of 
^(m) rpj^g higher the number of terms, the greater the 
possibility that spurious roots will appear. In practice, 
g ~ 30 is the highest order two-point Pade approximant 
that is computable for vortex solutions to the NLSE. 

In the next section, wc will demonstrate an alternative 
method that does not suffer from the limitations of the 



two-point Pade approximant. However, the Pade approx- 
imant is worth retaining because it provides interpolating 
functions in the form of rational polynomials which can 
reproduce the small and large x behavior of the wave- 
function to very high order. We note that a thorough 
treatment of the use of Pade approximants in the study 
of vortex and other solutions to the NLSE has been made 
byN. G. Berloff ,651]. 



3. Solution by Numerical Shooting 

The 2D NLSE in the form given by Eq. Qcan be 
solved by shooting. In this standard method jBa], one 
chooses the values of ?7(xo) and rj'{xo) for Xo 1. In 
this way, one can obtain an accurate relation between 
'?(xo) and 77'(xo) via the power series of Eq. (|A1|) . One 
then integrates the NLSE by initial value methods to- 
wards arbitrarily large values of x- The correct initial 
value of Tj and -q' leads to the vortex solution. Because 
the vortex solution lies on the boundary between diver- 
gent and nondivergent behavior, it is quite easy to tell 
when one has made a wrong choice: either 77 diverges to 
infinity or it oscillates and approaches zero. One chooses 
an initial value of ag, then iterates. Note that the bound- 
ary cannot be chosen at xo = 0, since 77 = is a valid 
solution to the 2D NLSE. Rather, a value of xo which is 
exponentially small is used, so as to ensure the accuracy 
of the power series solution. 

In practice, the order of the power series is never a 
limiting factor. For example, we worked with 40 terms. 
Up to ten additional terms were tried without finding 
any difference in the results. Instead, the two adjustable 
parameters in the calculation were xo and the number of 
digits of internal precision used in our numerical routine. 
Of these, it was the latter that most strongly affected 
the critical value of cq. In order to determine to the 
highest possible precision, it was necessary to use num- 
bers of higher than double precision. It was found that 
35 digits was a practical maximum for our computing 
capabilities. In each case, the number of digits of preci- 
sion was determined by comparing the results using 32, 
33, 34, and 35 digits of internal computational precision. 
In Table ^1 are shown the results. Note that they are 
consistent with and greatly improve upon those of the 
two-point Pade approximant shown in Table ^ As in 
Sec. lA 21 the calculations proved more computationally 
difficult at higher winding number. The reader may ask 
why obtaining such high precision in the value of Cq™'' is 
desirable. The reason is that each digit of precision brings 
the solution a few units of x closer to the exact vortex 
solution. In order to obtain a solution which is exactly 
on the boundary between divergence and non-divergence 
an infinite number of digits are required. We illustrated 
this extremely sensitive dependence of the determining 
coefficient on the number of digits of precision in Fig. ^ 
Figure ^ depicts three values of ao which approximate 
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FIG. 9: (color online) Approaching the vortex solution. The 
position of the first node xi is shown as a function of the num- 
ber of digits of precision k in the determining coefficient Qq™' 
for quantum vortex solutions to the 2D nonlinear Schrodinger 
equation, as defined in Eq. I)A14|I . Observe that all curves are 
close to parallel and linear in k, and as the winding num- 
ber m increases, the curves converge. Note that all azxes are 
dimensionless. 



ring soliton solutions, as discussed in Sec. II VI 

One finds an intriguing relationship between the posi- 
tion of the first node and the number of digits of precision 
k in the critical determining coefficient aQ™-*, where k is 
defined by Eq. (jA14p . In Fig. |^ is shown the position 
of the first node Xi(^)- One observes that the relation- 
ship is linear. For all values of the winding number but 
TO = 0, the curves lie nearly on top of each other, and all 
are parallel. Clearly, from Eq. Q, for large x the term 
which depends on to becomes negligible. Note that the 
use of the special case to = ensures that, in at least one 



case, the exact value of 



is known. The best values 



of Gq tor m ^ are given in Table Ull 

4. Modified Numerical Method for the Infinite 
Well 



In order to quantize the solutions in an infinite well 
in two dimensions, one holds the normalization and the 
cylinder radius to be constant. The form of the wave- 
function and the chemical potential can be obtained as 
follows. One calculates the dependence of the normaliza- 
tion on flo as 



(m) 



to 



k digits of precision, where k is defined by 



ao 



(1 - 10"''), 



(AM) 



In Eq. I)A14() . a small subtraction is made in the kth 
digit, so that a convergent rather than divergent solution 
is obtained. In Fig. QJa), the divergent solution is also 



depicted for the same value of fc, with ao = 



(1 



10~'°): i.e., a small addition is made in the fcth digit. It 
is in this sense that the vortex solution is a boundary 
between convergent and divergent solutions. 

The figure shows the solution obtained via numerical 
shooting for k = 16, 8, and 2 in panels (a), (b), and (c), 
respectively. A winding number of m = 1 is assumed. 
In (b) and (c), the usual Bessel function solution to the 
two-dimensional linear Schrodinger equation is shown for 
comparison. One sees that the higher the precision, the 
further the first node is pushed out towards large values 
of X- To move the node to infinity, an infinite number of 
digits of precision is required. All of the solutions except 
the divergent one depicted in Fig. ^a) are examples of 



Xj(ao) 



(A15) 



where Xji^-a) is the distance to the jth node in r](x)- 
As evident in Figs. ^ and |51 a more useful variable to 
determine the dependence of Afj on oq is the number of 
digits of precision k, as defined by Eq. HA14() . Note that k 
is not restricted to an integer value. The value of k also 
determines Xj- From Eq. H12|) . the chemical potential 
scaled to the energy associated with the cylinder radius 
is 



(A16) 



where fij is the chemical potential for the (j — l)th excited 
state, with j = I giving the ground state. The function 
liiij{Mj) can be calculated from Eqs. (|A15p and ljA16p . 
This gives the chemical potential as a function of the 
atomic interaction strength. 

Quantized modes for the infinite spherical well in three 
dimensions are calculated by similar methods. 
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